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Abstract. We give an analysis of a continuation algorithm for the numerical solution of the force- 
based quasicontinuum equations. The approximate solution of the force-based quasicontinuum 
equations is computed by an iterative method using an energy-based quasicontinuum approximation 
as the preconditioner. 

The analysis presented in this paper is used to determine an efficient strategy for the parameter 
step size and number of iterations at each parameter value to achieve a solution to a required 
tolerance. We present computational results for the deformation of a Lennard- Jones chain under 
tension to demonstrate the necessity of carefully applying continuation to ensure that the computed 
solution remains in the domain of convergence of the iterative method as the parameter is increased. 
These results exhibit fracture before the actual load limit if the parameter step size is too large. 



1. Introduction 

Quasicontinuum (QC) approximations reduce the computational complexity of a material simu- 
lation by reducing the degrees of freedom used to describe a configuration of atoms and by giving 
approximate equilibrium equations on the reduced degrees of freedom [5-7,9-11,13-16,18,20,21, 
23-25]. For crystalline materials, there are typically a few small regions with highly non-uniform 
structure caused by defects in the material which are surrounded by large regions where the local 
environment of atoms varies slowly. The idea of QC is to replace these slowly varying regions with 
a continuum model and couple it directly to the atomistic model surrounding the defects. The ma- 
terial's position is described by a set of representative atoms that are in one-to-one correspondence 
with the lattice atoms in the atomistic regions but reduce the degrees of freedom in the continuum 
regions. 

Quasi-static computations in material simulations explore mechanical response under slow ex- 
ternal loading by fully relaxing the material at each step of a parameterized path of external 
conditions. Such simulations can model nano-indentation, stress- induced phase transformations, 
and many other material processes. The characteristic feature when using this technique is that 
the process to be modeled occurs slowly enough that dynamics are assumed to play no role in 
determining the relaxed state. This paper focuses on applying continuation techniques [4,8,12] to 
the nonlinear equilibrium equations of the force-based quasicontinuum approximation (QCF). 

1.1. Choosing a Quasicontinuum Approximation. There are many choices available for the 
interaction among the representative atoms, especially between those in the atomistic and contin- 
uum regions, which has led to the development of a variety of quasicontinuum approximations. Cri- 
teria for determining a good choice of approximation for a given problem are still being developed. 
Algorithmic simplicity and efficiency is certainly important for implementation and application, but 
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concerns about accuracy have led to the search for consistent schemes. Since the forces on all of 
the atoms in a uniformly strained lattice are zero, we define a quasicontinuum approximation to be 
consistent if there are no forces on the representative atoms for a lattice that has been deformed by 
a uniform strain. We note that the atomistic to continuum interface is typically where consistency 
fails for a QC approximation [7,9,18,24] as is the case for the original QC method [25]. 

For static problems, QCF is an attractive choice for quasicontinuum approximation because it is 
a consistent scheme [7, 18] that is algorithmically simple: the force on each representative atom is 
given by either an atomistic calculation or a continuum finite element calculation. The algorithmic 
simplicity of the force-based quasicontinuum method has allowed it to be implemented with adaptive 
mesh refinement and atomistic to continuum model selection algorithms [1-3,6,18,19,22]. The 
trade-off for the consistency and algorithmic simplicity of QCF is that it does not give a conservative 
force field, although it is close to a conservative field [7]. Thus, QCF is a method to approximate 
forces, rather than a method to approximate the energy. 

Quasicontinuum energies have been proposed that utilize special energies for atoms in an in- 
terfacial region [9,24], and corresponding conditions for consistency have been satisfied for planar 
interfaces [9]. However, there is currently no known consistent quasicontinuum energy that al- 
lows general nonplanar atomistic to continuum interfaces and mesh coarsening in the continuum 
region (other than the computationally intensive constrained quasicontinuum energy discussed in 
Section [2.2p . We will, however, investigate the use of quasicontinuum energies as preconditioners 
for the force-based quasicontinuum approximation. 

1.2. Solving Equilibrium Equations by Continuation. Our goal is to efficiently approximate 
the solution z(s) to the QCF equilibrium equations 

FQ^'^{z{s),s) + f{s) = 0, se[0,l], 

where z(s) € M" are the coordinates of the representative atoms that describe the material and 
the parameter s G [0, 1] represents the change in external loads such as an indenter position or 
an applied force. Using continuation, we start from z(0) which is usually easy to find (such as a 
resting position) and follow the solution branch by incrementing s and looking for a solution z(s) in 
a neighborhood of the previous solution. The continuation approach that we analyze in this paper 
has been used to obtain computational solutions to materials deformation problems [6, 18] and is 
implemented in the multidimensional QC code [17]. 

The approach that we will investigate in this paper uses constant extrapolation in s to obtain 
initial states for the iterative solution of F^^^ {z,, s) at a sequence of load steps Sq where = sq < 
■si < • • • < sq = 1 and solves the iterative equations using a preconditioner force F'^^^ that comes 

from a quasicontinuum energy f '^'^■^(z, s), that is, F^^^ (z, s) = —^^—q^^^- We will focus our 
analysis on a specific preconditioner later, but the splitting and subsequent analysis works for any 
choice of quasicontinuum energy. The outer iteration at a fixed step Sq is given by 

F«^^(zf 1, = F«^^(Z^, Sq) - F^izP, Sq) 

= -/(s,)-F^(zP,S,), p = 0,...,Pq-l, (1.1) 

where 

F^{z,s) ■.= f{s)+FQ^^{z,s) 

is the residual force and 

F^{z,s) :=F«^^(z,s)-FQ^^(z,.) 
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is sometimes called a "ghost force correction" in the mechanics literature [18]. We will consider 
preconditioner forces F^'-^^[z, s) that differ from F^'~"^{z, s) only in atomistic to continuum inter- 
facial regions so that the ghost force correction is inexpensive to compute. Since the preconditioner 
forces come from an energy, the outer iteration equations p.ip for Zq'^^ can be solved by an inner 
iteration that finds the local minimum in a neighborhood of Zq for 

«^^(Z, - {f{Sq) + F^iz^q, Sq)) ■ z] (1.2) 

using an energy minimization method starting with initial guess z = Zq. 

We give an analysis to optimize the computational efficiency of the continuation algorithm (II. ip 
by varying the parameter step size, hq = Sq — Sq-i, and the number of outer iterations, Pq. Our 
analysis first considers the goal of computing an approximation of z(s) uniformly for s G [0, 1] to a 
given tolerance, e. The proposed strategy selects the step size hq so that the initial iterates z° are 
in the domain of convergence of the outer iteration and so that the tolerance is achieved by the 
continuous, piecewise linear interpolant of the solution at each parameter Sq. The required accuracy 
at the Sq is achieved by the fast convergence of the iteration (jl.ip . As e — > 0, our analysis gives 
that hq ^ and — > oo for all q = 1, . . . ,Q, so that an efficient way to achieve increased accuracy 
uses a balance between small step size for accurate interpolation and a large number of iterations 
per step. 

We then consider the goal of efficiently computing the final state z(l) to a given tolerance, e. 
For this goal, the result of our analysis states that an efficient strategy fixes the number of outer 
iterations to Pg = 1 at all but the final step and takes the largest possible steps, hq, such that the 
initial guesses, z^, remain in the domain of convergence of the iteration. This strategy determines 
the number of steps, Q, independently of e. The required tolerance, e, is then achieved at s = 1 by 
doing sufficiently many iterations, Pq > 1. In this case, only Pq ^ co as e — > 0. 

We give numerical results for the deformation of a Lennard-Jones chain under tension that 
demonstrate the importance of selecting the step size, hq, and number of outer iterations, Pq, so 
that the iterates Zq remain in the domain of convergence of the iteration. The numerical experiment 
shows that our algorithm diverges (the chain spuriously undergoes fracture) if we attempt to solve 
for the deformation corresponding to a load near the limit load by a step size hi = 1. 

We give a derivation of the force-based quasicontinuum approximation and the energy-based 
preconditioner in Section [2j In Section [3l we analyze the equilibrium equations and their iterative 
solution. In Section [H we apply Theorem 13.11 to a Lennard-Jones chain under tension to obtain 
bounds on the initial state that guarantee convergence of our iterative method to the equilibrium 
state, to obtain convergence results for our iterative method, and to demonstrate the need for 
continuation by the showing that the chain can undergo fracture if we begin our iteration outside 
the prescribed neighborhood. Section [5] presents the continuation method and Sections [6] and [7] give 
an analysis to guide the development of an efficient algorithm using the quasicontinuum iteration. 
We collect the proofs of several lemmas in a concluding Appendix El 

2. Quasicontinuum Approximations 

This section describes a model for a one-dimensional chain of atoms and a sequence of approx- 
imations that lead to the force-based and energy-based quasicontinuum approximations. While a 
one-dimensional model is very limited in the type of defects it can exhibit, its study illustrates 
many of the theoretical and computational issues of QC approximations. 

We treat the case where atomistic interactions are governed by a pairwise classical potential (j){r), 
where (p is defined for all r > 0. A short-range cutoff for the potential is a good approximation 
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for many crystals, and in the following analysis we use a second- neighbor (next-nearest neighbor) 
cutoff, as this gives the simplest case in which the atomistic and continuum models are distinct [7]. 

2.1. The Fully Atomistic Model. Denote the positions of the atoms in a linear chain by yi for 
i = —M, . . . , M + 1 where yi < yi+i and denote the position where the right-hand end of the chain 
is fixed by yM+i = yM+i{s) for a parameter s € [0, 1]. The second- neighbor energy for the chain is 
then given by 

M M-l 

£"{y, s) := -y^)+ (A(y^+2 - Vi) (2.1) 

i=-M i=-M 

where y := {y-M, ■ ■ ■ ,yM+i)- We also assume that the chain is subject to an external potential 
energy which we assume for simplicity to have the form 

M 
i=-M 

Section U] describes a numerical example with a boundary dead-load given by /_a/(s) ^ and 
fi{s) = for all interior atoms i = —M + 1, . . . , M. 
We want to find local minima of the total energy, 

£Lar-=s''iy,s) + £Uy,s) (2.2) 

subject to the boundary constraint yM+i = yA.i+i{s). The equilibrium equation for the fully atom- 
istic system (j'2.'2p is given by 

F,^y{s),s) + fiis) = 0, i = -M,...,M, 

where the atomistic force is given by 

F-{y, s) := = [<A'(yi+i - y^) + <A'(y.+2 - y^)] 

- [(t>'{yi - Vi-i) + (P'ivi - yi-2)] , 

for i = —M, . . . ,M where the terms </>'(yi — yj) above and in the following should be understood 
to be zero for i ^ {— M, . . . , M + 1} or j ^ {— M, . . . , M + 1}. In the remainder of this section, we 
will not explicitly denote the dependence on the parameter s. 

2.2. The Constrained Quasicontinuum Approximation. The constrained quasicontinuum 
approximation finds approximate minimum energy configurations of (j2.2p by selecting a subset 
of the atoms to act as representative atoms and interpolating the remaining atom positions via 
piecewise linear interpolants in the reference configuration. We denote by ao the ground state 
lattice constant for the potential (j){r) with a second-neighbor cutoff, that is, 

ao := argmin</>(r) + (j){2r) 

(see Section [3] or [7]). We then set the reference positions of the atoms in the chain to be 

Xi := iao for i = —M, . . . , M + 1. 

We let Zj := yi^ denote the representative atom positions where j = —N, . . . , -|- 1, and where 
i-N = —M, In+i = M + 1, and ij < ij+i- We are interested in developing methods for <C M. 
We can obtain the positions of all atoms yi from the positions of representative atoms zj by 

N+l 

yi{z) = Y Sj{x^)zj for i = -M, . . . , M + 1, 

j=-N 
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where z := {z-n, ■ ■ ■ , -^Af+i) and the Sj{x) are the continuous, piecewise hnear "shape" functions for 
the mesh constructed from the reference coordinates xi . of the representative atoms, more precisely, 



Sj{x) :-- 



'O, ifx<X£^._j, 

{x - xe^_^)/{xi^ - xi^_^), if xi^_^ < X < xe., 

- x)/{xi^^^ - xe.), if xe^ < x < Xi^^-,, 

0, iix>xe^^,. 

The constrained quasicontinuum energy is then given by 

f^«^(z) := r(y(z)), 

and the constrained external potential energy is given by 

£:£?^(z) :=^,",,(y(z)). 



(2.3) 



Using (12. 3p and the chain rule, we obtain the conjugate atomistic force, that is, the force on the 
reduced degrees of freedom induced by the atomistic forces. We find that 



dzj 



for j = —N, . . . ,N, and the conjugate external force is given by 

j=o '■'^^ ' 1=1 

for j = —N, . . . ,N, where 

Uj := ij+i - ij 

is the number of atoms between zj and Zj+i (the end atoms are only counted half). The equilibrium 
equations for the total constrained quasicontinuum energy, 

^^?a?(-) :=£:^«^(z)+^fi^(z) 

are then given by 

«^(z) + = 0, j = -N,...,N. 

The constrained quasicontinuum approximation is attractive since it gives conservative forces 
and since it is the only known quasicontinuum energy that is consistent when generalized to multi- 
dimensional approximations [9]. The constrained quasicontinuum approximation is also attractive 
since its conjugate forces (|2.4p are located at only 2N representative atoms; however, we must still 
compute the forces at all 2M atoms which makes it computationally infeasible. Some computational 
savings can be made in the interior of large elements by separating the energy computations into 
element energy plus surface energy; however, in higher dimensions the large number of atoms near 
element boundaries makes the constrained quasicontinuum approximation impractical. Finally, it 
is attractive because its approximation error comes only from the restriction to linear deformations 
within the element making it possible to analyze using classical finite element error analysis. 
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2.3. The Local Quasicontinuum Energy. We now recast the constrained approximation in 
terms of continuum mechanics to introduce the local quasicontinuum energy which is simply a 
continuous, piecewise linear approximation of a hyperelastic continuum model where the strain- 
energy density is derived from the atomistic potential, (?!>(r). This energy efficiently approximates 
the conjugate force at the representative atoms. We have that [7] 

TV 

£CQC{z)= ^ LjW{Dj)+Sb{D_N) 

j=-N 

TV 

+ Yl S{Dj_i,Dj)+Sb{DN), 

j=-N+l 

where 

L,:=.,^^,-.,^ a^d Dr- " ^ 



xe,+, - XI. 

are the length and deformation gradient of the jth clement, and 

^^^^ ^ (t>{Dao) + ct>{2Dao) ^ 

is the strain-energy density for an infinite atomistic chain with the uniform lattice spacing DaQ. 
Here [7] 

Sb{D) = - ^0(2Dao), 

S{Di,D2) = - ^(f>{2Diao) + 0(L'iao + £'200) - ^^{2D2ao), 

can be considered to be a surface energy and an interfacial energy respectively. 

Since S{Dj-i, Dj) is a second divided difference, the interfacial energy is small in regions where 
the strain is slowly varying. We obtain the local quasicontinuum energy by neglecting the interfacial 
energy and surface energy to obtain 

N 
j=-N 

and we have the corresponding conjugate atomistic force 

^'«-^-^<-^)-^(--). 

We note that -^/'(z) depends only on Zj^i, Zj, and Zj^i. This approximation is computationally 
feasible since the work to compute all the forces is proportional to N. The approximation error now 
has two components: the linearization within each element that is inherited from the constrained 
quasicontinuum approximation plus the operator error incurred by ignoring interfacial terms. In 
cases where the deformation gradient Dj is slowly varying on the scale of the representative atom 
mesh, both sources of error will be small and the local approximation will be highly accurate, as 
is expected for a sufficiently refined finite element continuum model. Mesh refinement can be used 
to reduce both sources of error, but even mesh refinement to the atomistic scale cannot remove 
the interfacial error in the neighborhood of defects since the deformation gradient varies rapidly on 
the atomistic scale. Thus, the atomistic model must be retained near defects to obtain sufficient 
accuracy. 
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2.4. The Force-based Quasicontinuum Approximation. We can obtain a quasicontinuum 
approximation that is accurate in regions where the deformation gradient Dj is rapidly varying, 
such as in the neighborhood of defects, and maintains the efficiency of the local quasicontinuum 
method by combining them in the force-based quasicontinuum approximation (QCF). In QCF, we 
partition the chain into "atomistic" and "continuum" representative atoms and define the force on 
each representative atom to be the force that would result if the whole approximation was of its 
respective type, that is. 



?QCF, 



z) : = 



f;(z) 

if (z) 



if representative atom j is atomistic, 
if representative atom j is continuum. 



(2.5) 



With this convention, for example, the forces on a continuum representative atom are determined 
solely by the adjacent degrees of freedom regardless of how close any atomistic representative atoms 
may be. 



---Q 



O 



ZK-1 



ZK 



ZK+1 



ZK+2 



ZN + 1 



Figure 1. One end of the quasicontinuum chain, highlighting the interface. Filled 
circles are atomistic representative atoms, whereas the unfilled circles are continuum 
representative atoms. 

For simplicity, we will consider a single atomistic region symmetrically surrounded by continuum 
regions large enough that no atomistic degrees of freedom interact with the surface. We let the rep- 
resentative atoms in the range j = —K + 1, . . . , X be atomistic and in the ranges j = —N, . . . , —K 
and i^T + 1, . . . , A'^ be continuum. Figure [1] depicts one end of the quasicontinuum chain. We note 
that the atomistic model has surface effects at the two ends of the chain, but the local quasicon- 
tinuum model does not have surface effects. Thus, this arrangement of representative atoms with 
continuum representative atoms at the ends of the chain will not give surface effects within the QC 
approximation . 

We assume that Vj = 1 for j = —K — l,...,K + \. This guarantees that Vj = 1 within the 
second-neighbor cutoff radius of any atomistic representative atom and thus allows -F"^^^(z) and 
F^(z) to be computed without interpolation. The forces are then given by 

+ 2^'{2r,)] 

/(r,-_i) + 2</>'(2r,_i)], 

+ 0'(r,-+r,+i)] (2.6) 



(z) 



-[<? 



-N + l<j< -K, 



'(rj-_i)+</)'(rj_ 
+ 20'(2r,)] 



,i+rj_2)], -K + l<j<K, 



where 



Djao 



_i) + 20'(2r,_i)] 

{Zj + i - Zj) 



K + 1 < j < N, 



-N, 
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is the deformed lattice spacing within the jth. element. 

2.5. An Energy-Based Quasicontinuum Approximation. There are many quasicontinuum 
energies [9,18,24] that can be used to precondition the iterative solution of the force-based quasi- 
continuum approximation (jl.ip . We will give an analysis and numerical experiments for the quasi- 
continuum energy described in [18] and denoted here by S^^^ because it seems to be the simplest 
to implement and because it converges sufficiently rapidly. Here and in the following QCE will refer 
specifically to the energy described in [18] , whereas in the introduction it represented any possible 
choice of quasicontinuum energy. 

QCE assigns an energy to each degree of freedom according to the model type (atomistic or 
continuum), and the sum of all such energies gives the total QC energy for the chain. We use the 
same distribution of atomistic and continuum representative atoms as above. Then the atomistic 
representative atoms, located in the range j = —K -|- 1, . . . , have energy given by 

£^{z) := ^ [0(r^-) + 0(rj- + rj+i) + 0(rj-i) + 0(rj_i + r,„2)] , (2.7) 

and the continuum representative atoms, located in the range j = —N,...,—K and j = K + 
1, . . . , -|- 1, have energy [7] given by 

£^{z) := i [L.WiD,] + L,_iVF(I),_i)] (2.8) 

where the energy density W{Dj) is considered to be zero for j < —N or j > N. The quasicontinuum 
energy, £'^^^{z,), for the chain is given by 

-K K N+l 

fQ^^(z) = £:/(z) + (^) + E ^/(^)- (2.9) 

j=-N j=-K+l j=K+l 

In (|2.5p . we assign forces according to representative atom type whereas here we have assigned a 
partitioned energy according to representative atom type. 

We now mention other QC energies, although they will not be used in the following. In [24], 
the quasinonlocal method is proposed to attempt to remove the interface inconsistency by defining 
a new QC energy. For this method, special interface atoms are defined that behave in a hybrid 
fashion, interacting atomistically with a neighbor if that neighbor is atomistic, but using the local 
approximation to determine the interaction energy otherwise. For example, if we denote represen- 
tative atoms j = K and X -|- 1 to be quasinonlocal, then their energy would be 

Sf{z) = ^ [^{rj) + ^{2rj) + 0(r,_i) + 0(r,_i + r,-_2)] . (2.10) 

However, this method only gives a consistent quasicontinuum energy for a limited range of in- 
teractions (second- neighbor in one dimension), and further inconsistencies are introduced when 
attempting to coarsen the continuum region in higher dimensions [9]. A more general approach 
that applies to longer-range interactions is given in [9], but this approach to the development of 
consistent quasicontinuum energies is also currently restricted to planar interfaces in higher dimen- 
sions. 

3. Convergence of the Iterative Method to Solve the QCF Equations 

We now give a theorem for the convergence of the iterative algorithm (jl.ip to solve the QCF 
equilibrium equations. Specifically, we give a domain in which the iteration is well-defined and 
a contraction. In the following, this contraction will be an essential portion of the continuation 
method that is applied to solve the final equilibrium problem. 
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Our result extends the theorem in [7] by allowing the removal of the hypotheses on the external 
force, f := (/-m, • • • , /m), by utilizing mixed boundary conditions in the problem analyzed in this 
paper rather than the free boundary conditions analyzed in [7] (the different assumptions lead to 
different constants in the inequalities). In this section, the dependence on s of both the solution, 
r, and external force, f , is again suppressed. 

Since the QCF forces ()2.6p and the QCE energy (j2.9p depend only on the interatomic spacings, 
{rj}jL_j^, the analysis of the iteration is simplified by formulating the problem in terms of forces 

on the lattice spacing, {rj}jL_j^, rather than on representative atom positions, {zj}^J'_}j^. We note 
that since ^^r+i = ijM+iis) is fixed, there is a one-to-one mapping z ^ r. For the energy-based 
quasicontinuum approximation, we define ^^"^^(r) to be the force conjugate to the representative 
atom spacing Zj+i — z 



3 

j, namely {r 



-u- ^ ^^2^, (z). This conjugate force satisfies 



-N <j < -K 
j = -K- 1, 



+ i cP' (rj + r, +1 ) + 0' (2r,- ) , j = -K, 

(p'irj) + ^(j)'{rj + rj-_i) + (P'{rj + r^+i), j = -K + 1, 
<i6'(r,) + (/)'(r, + r,_i) + 0'(r,- + r,+i), -K + 2<j <K-2, 



We have from the chain rule that [7] 



dzj 
I QCE, 



(^) 



^wo.(r)+V;«5;^(r), j = -iV, . . . , iV, 



where := 0, so it follows by summing (j3.ip that 



(3.1) 



i=-N 



If we define an analogous quantity ijj^^^ (r) by setting 



-N. 



i=-N 



then we have that 
-^?^^(r) 



.N. 



, . . . , T , 



-N <i < -K, 
-K + 1 < j < K, 
K + 1 < j < N, 



(3.2) 



0'(r,) + 20'(2r,-), 

(p'irj) + (l)'{rj + rj_i) + (l)'{rj + r^+i) + I.k, 

[cf>'{r,) + 2<P'{2rj)+I_K-lK, 

where Ij = 2(p'{2rj) — (j)'{rj + rj_i) — (j)'{rj + ^'j+i). The external force is likewise made conjugate 
to the representative atom spacing by summing 

3 



E/^ 



3 



(3.3) 



-N 
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It follows from (j3.2p and (|3.3p that a configuration z is a solution to 

F^^^{z)+fj = 0, j = -N,...,N, (3.4) 
if and only if the corresponding r is a solution to 

V^f ^^(r) + CD, = 0, j = -N,...,N. 

We will iteratively solve the equilibrium equations (j3.4p by using F^^^ as a preconditioner for 
pQCF_ The convergence theorem we prove later in this section means that QCE is quite close to 
QCF, and the iterative equations converge rapidly. Thus, if we use standard energy minimization 
algorithms to solve the QCE equations at each iterative step and utilize the fast convergence of 
the QCE solution to QCF, then we get an efficient solution method for the QCF equations with its 
inherent advantages of consistency and simplicity. The iterative equations are 

+ V^f (r^') + CD, =0, j = -TV, . . . , iV, (3.5) 

where the correction force is 

V^,«(r):=V',«^^(r)-V',«^^(r). (3.6) 

3.1. Assumptions on the Atomistic Potential, (/>(r). Before stating our result about the con- 
vergence of the iteration (j3.5p . we make explicit the assumptions on the potential, A prototypical 
function fitting these assumptions is the Lennard-Jones potential, 

<t^i^) = jr2 - y,- (3-7) 

We recall that the energy density corresponding to (j){r) for the second-neighbor energy (j2.ip is 
given by W{D) = ^ {(l){Dao) + (/)(2Z)oo)) where ao is the equilibrium bond length of a uniform 
chain, that is, it is the minimum of (j){r) + (f)[2r). 

We will assume that (j){r) G ((0,00)) and that it satisfies the following properties that are 
illustrated in the Lennard-Jones (j3.7p case in Figures [2] and [3l There exist fi, f2, and D such that 

(f)" {r) > for < r < fi and (j)"{r) < for r > fi, 
(/)"'(r) < for < r < f2 and (/>"'(r) > for r > f2, 
W'{D) < for < L> < 1 and W'{D) > for D > I, 

W"{D) > for < D < 5 and W"{D) < for D > D, 

< ao < fi < f2 < 2ao, 

1 < D. 

We note that D is the deformation gradient of a uniform chain at the load limit. 

The following theorem gives sufficient conditions on the existence of a region r = (r_Ar, . . . , r/vr) G 
i7 = {rL,ru)'^^~^^ in which the iteration (13. 5p is well-defined and a contraction. We see that 
under these conditions QCE is an efficient preconditioner for the force-based equations, giving a 
contraction mapping for the iteration. The idea is that these quasicontinuum approximations are 
quite close, so that the solution of QCE gives a good approximation to the solution of QCF. 
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Figure 2. The Lennard- Jones potential (|3.7p demonstrates the prototypical behav- 
ior of ^(r) and its derivatives. 



w w w 




Id Id Id 

Figure 3. The energy density, W{D), corresponding to the Lennard- Jones poten- 
tial (j3.7p and its derivatives. 



Theorem 3.1. For a given conjugate external force, {^-m^ ■ ■ ■ ,^n) suppose that there exist r^ 
and rjj such that 

y < r-L < ru, (3.8) 
0"(rt/) + 21(/)"(2ri) > 0, (3.9) 
(t)\rL) + 60'(2rL) - 40'(2rc/) < < ^'{ru) + 64>'{2ru) - 40'(2ri), (3.10) 
for j = —N, . . . ,N. Then for every G := (r^, r^/)^^^^ there is a unique r^+^ £ 0, such that 

^QCE^^P+i^ + (r^) + =0, j = -N,...,N. (3.11) 
We also have that the induced mapping is a contraction: if — > r^+^ and — > s^"*'^, then 



I oo 



<P"iru)-5\(P"{2rL)\ 



loo ' 



where we have from (j3.9p t/iai 

16|(/."(2ri)| 



< 1. 



0"(r[/)-5|0"(2rL)I 

We start by remarking on the theorem's assumptions. The second inequality in (j3.8p states 
that is acting as a lower bound on miuj rj and ru as an upper bound. The first inequality is 
chosen for convenience so that (j)"{2r) is monotone. (We note that for Lennard- Jones and similar 
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potentials, it is physically a very reasonable assumption due to the stiffness of the interactions.) 
Condition (j3.9p ensures diagonal dominance of the Jacobian matrix for F^'-^^ and ensures that the 
contraction estimate is less than 1. Finally, the condition on in (j3.10p restricts the external 
forces sufficiently to allow a simple degree theory argument to prove existence of solutions. 

It is possible to choose a fairly large range for r when the external forces are far from the load 
limit of the chain. However, as maxj approaches the load limit, r must approach the tensile limit 
which makes the estimates much more sensitive. This reduces the size of {rL,ru)- The hypotheses 
of Theorem 13.11 guarantee that the iteration (j3.1ip converges to the QCF approximation of a stable 
atomistic solution. 

This theorem is a modification of Theorem 5.1 in [7]. The proof there models the QCE equations 
as a perturbation of the fully local quasicontinuum energy, S^. The proof here follows by modifying 
the original proof to handle the new terms that arise from removing the assumption of symmetry 
of . These new terms can be estimated by the techniques used to estimate similar terms analyzed 
in [7]. 



4. The Deformation of a Lennard-Jones Chain under Tension 

In this section, we consider the application of Theorem 13.11 to a chain modeled by the Lennard- 
Jones potential ()3.7p . The deformation of the fixed end, y^j^i(s), can be set arbitrarily since the 
dependence on yM+iis) is given by a uniform translation. To obtain a uniform tension, we model 
the external force for the fully atomistic chain by /-a/ = — ^ and fj = for j = —M + 1, . . . , M. 
It then follows from (13. 3p that the conjugate external force for the QC approximation is given by 
$j = $ for ah j = -N, ...,N. 

There are uniform solutions to the QCF equations up to the load limit ^max = 2.7810, that is, 
if (p'{r) + 2(p'{2r) = then *0'^^(re) = ^>e, where e = (1, 1, . . . , 1) G M^^+i and *Q'^^(r) := 

(i^-N^ • • • ) V'^*^^(i")^ • We apply an external force very close to the load limit to get an example 
where continuation is necessary to ensure that the preconditioned equations converge. Define the 
loading path $(s) = 2.76,se. Then solutions r(s) to the QCF equations satisfy r(s) = r(s)e where 

0'(r(s)) + 2(p'{2r{s)) = 2.76s, s G [0, 1]. (4.1) 

For any s G [0, 1], we can apply Theorem 13.11 to this example by picking and rjj such that 

r{ru)-5\<P"{2rL)\ ' ^ ' ' 

to conclude that the iterative equation using the QCE preconditioner is a contraction mapping 
with contraction constant a, provided that (|3.8p - (|3.10p holds. We find ri and ru symmetrically 
positioned about r{s) by substituting ri(s) = r(s) — 6{s) and ru{s) = r(s) + 6{s) into (|4.2p with 
r(s) given by the solution to (|4.ip to obtain 

cj)"{r{s) + 6{s)) + (5 + 16/a)(l)"{2{r{s) - S{s))) = 0. (4.3) 

It can be checked that (j3.8p - (j3.10p are satisfied with <I>j = $ = 2.76s for j = —N, . . . ,N. Therefore, 
for any initial guess r'^ G [r(s) — 6{s),r{s) + 5{s)]'^^~^^, the iteration step (13. Sh is a contraction 
mapping for all n with contraction rate a and limit point r(s) = r(s)e. Figure [4] depicts the solution 
r(s) along with four contraction intervals that correspond to a = ^, ^, ^, and |. For every a < 1 the 
corresponding 6{s) is decreasing. In Section [7] we consider the contraction region corresponding to 
a = i since this contraction region terminates just beyond our maximum applied load, (5(1.001) = 0. 
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s > s > 



(a) (b) 

Figure 4. (a) Loading response, r{s), for the Lennard- Jones chain is surrounded 
by contraction regions, r(s) it S{s), corresponding to a = |, |, ^, |. The contraction 
constant a increases with distance from r(s). (b) Detail shows the contraction region 
in a neighborhood of s = 1. 

4.1. Fracture at the Interface. To demonstrate the need for continuation methods for the above 
example, we describe the performance of a modified version of qcld, a code by Ellad Tadmor for 
solving the QCF equations using QCE as a preconditioner with a nonlinear conjugate gradient 
method for solving the inner iteration. We attempt to directly solve (jl.2p starting from the energy- 
minimizing lattice spacing and using only a single loading step, Q = 1, which gives the following 
minimization problem. We have 

r^ = argmin [*<3Ci?(j.) _ (2.76 + *'^(ao)) • R(r)] , (4.4) 

r 

where R(r) = (z^_Arr_Ar, . . . , i^n^n) denotes the representative atom spacing. 
We consider an uncoarsened QC chain, with 

M = N = 7, 

undergoing external loading as described in Section HI The chain is partitioned with 

K = 3 

which means that there are six atomistic representative atoms surrounded symmetrically by two 
groups of five continuum representative atoms. The QCF solution r(s) given by (14. Ih and the 
contraction parameters 5{s) and a given by ()4.3p do not depend on the size of the chain; however, 
the QCE preconditioner solution will depend on the size and composition of the chain because it has 
a non-uniform solution due to the atomistic to continuum interface. Because the exact solution is 
a uniformly deformed chain, our problem is unchanged by any coarsening of the continuum region. 
While this does not illustrate the power of QC approximations to reduce computational complexity, 
it provides a simple case in which to analyze loading up to a singular solution, in this case fracture. 

The chain fractures in the atomistic to continuum interface. In the interface, QCE behaves like 
a continuum material with varying stiffness which is why it fails to be a consistent scheme. The 
corrections ^'^(ao) act to counterbalance this effect by compressing the high stiffness regions, and 
adding tension to the low stiffness regions. Fracture occurs due to the fact that the corrections 
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Figure 5. A close-up of the atomistic to continuum interface showing fracture 
that occurs when continuation is not employed. The three layers represent three 
steps of a single conjugate gradient minimization for the iterative equations (j4.4p . 
where the position of zq has been normalized to align the chains. The upper layer is 
the initial, undeformed state r = ag. The middle layer shows a subsequent iteration 
where the chain is nearly uniformly deformed and close to the QCF solution. The 
lower layer shows later iteration where a clear separation of atom pairs occurs. None 
of the states shown is a solution to the minimization problem, and the numerical 
algorithm eventually terminates without finding a minimum. 



()3.6p applied in the atomistic to continuum transition are a model correction at the equilibrium 
bond length, oq, but are much too strong at the stretched configuration. The overcorrections add 
to the very large external force and exceed the load limit for the QCE chain (see Figure [5]). The 
above estimates show that the continuation method described in Section [5] provides a convergent 
method for computing the deformation of the chain at the load <I> = 2.76 with the qcld code. 

5. Solution of the QCF Equations by Continuation 

In this section and the following two, we give an analysis of the solution of the QCF equilibrium 
equations by continuation. We will present our results in a general setting that focuses on the 
contraction property of the preconditioned equations. Because we only use the abstract contraction 
property, the continuation analysis given here will apply to higher dimensional QC approximations 
provided one has a contraction result similar to Theorem 13.11 Given G € C^(M"'~'~^; M"), our goal 
will be to approximate a curve of solutions r(s) G (^^([0, to 

G(r(s),s) = for sG [0,1], (5.1) 

where 

det VrG(r(s),s) / for sG [0,1]. 
We will later apply this theory to QCF by considering the solution of the equations 
G(r(s), s) := i;'^^^{r(s),s) + $(s) = for all s G [0, 1]. 
Let k{s) be a bound on r'(s), that is, ||r'(s)||QQ < k{s), which gives 

||r(t) - r(s)||^ < / /c(r)dT for alH > s. 

J s 
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We assume further that for each s G [0, 1] there is an iterative solver Tg : M" that is locally a 

contraction mapping with fixed point r(s). That is, there is an a < 1 such that for every s G [0, 1], 
there is a radius 6{s) with the property 

- pIIoo ' - qlloo < ^is) \\TsP - T^qlL < « llp - qlloo • 

We saw in Section[3]that such a radius 6{s) can be obtained for Ts given by the iterative method ()3.1ip 
if the hypotheses of Theorem 13.11 are satisfied. 

Let = sq < si < • • • < sq = 1 be a sequence of load steps where at each point Sq we wish to 
compute Tq, an approximation to r{sq). Beginning from an initial guess r^, the iterative solver T<j^ 
is applied to (jS.ip . keeping Sq fixed. This generates a sequence of approximations = Tf^r^ for 
p = 1, . . . , Pq, where Tf^ denotes p compositions of the operator Tg^ and Pq denotes the number of 
iterations at step m. We then let = rg'. The choice of initial guess rjj is typically made using 
polynomial extrapolation, and here we choose = rq_i, which is zeroth-order extrapolation. 

We will now give an analysis of how to choose the load steps = sq < si < ■ ■ ■ < sq = 1 and the 
corresponding number of iterations Pi, . . . ,Pq to efficiently approximate r(s) with respect to two 
different goals. We first consider the efficient approximation of r(s) in the maximum norm for all 
s G [0, 1], and we then consider the efficient approximation of the end point r(l). We note that our 
analysis only gives an upper bound for the amount of work needed to compute an approximation of 
our chosen goal to a specified tolerance since we use a uniform estimate for the rate of convergence 
a rather than the decreasing a as we converge to the solution that we can obtain from Theorem 13. II 
and displayed in Figured! 



6. Efficient Approximation of the Solution Path in the Maximum Norm 

For simplicity, we will first consider a uniform region of contraction radius 6{s) = 5, a uniform 
bound k{s) = k, a uniform step size h = 1/Q = Sq — Sq-i, and a uniform number of iterations at each 
step Pq = P. We will denote the continuous, piecewise linear interpolant of r(sg) G M", q = 0, . . . ,Q, 
by 2r{s); and we will denote the continuous, piecewise linear interpolant of G M", q = 0, . . . ,Q, 
by f(s). We will determine an efficient choice of h and P to guarantee that 

max ||r(s)-f(s)||^ <2e, (6.1) 

sG[0,l] 

where we assume for convenience that 2e < 6. 

We will assume that r(s) G C^([0, 1];M"'), so there exists a constant k2 > such that 

max ||r(s) — Xr(s)|l < koh"^. 
SG[0,1] °^ 

We can then ensure that 

max ||r(s) -Ir{s)\\^ < e 
se|o,i] 

by choosing h < \J ejki- We can thus satisfy ()6.ip by guaranteeing that 

max ||r(sg) -rgll^ < e. (6.2) 

Now if ||r(sq„i) - rg_i||^ < e, then 

- A\ loo - ll^^^?) " r(«'/-i)lloo + I|r(s9-i) - rq-ill^ 
< A:/i + e. 
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We choose < h < so that is in the region of contraction 

||r(s„)-r°|| <kh + e<5. 

M ^ y Moo — — 

We then choose P to achieve the desired error 

< a^{kh + e) < e. 

We can thus guarantee that ||r(sq) — r^||^ < e by doing P iterations where 



In 

m = - 



i+kh 



In a 

The computational work to obtain (j6.2p can then be bounded by 



W{h) = —-^ = ^ for < /i < 

h Aim a 



We have by the Mean Value Theorem that 



dWih) k 



dh hlna 

for < /i< 



ln{e + kh)- In e 1 



kh {e + kh) 



< 



k ■ _ 

We can finally obtain (|6.1|) by choosing 



hopt = min | \/e7^ 



min|ln|,ln J^^' A 

P{hopt) = ] ^ oo as e ^ 0. 

ma 

As e goes to zero, the second criterion becomes active so that the step size is determined by the 
interpolation estimates rather than the size of the contraction region. The number of steps grows 
as and the number of iterations per step grows like like 2 hi a • 

7. Efficient Approximation of the Solution at the Final State 
In this section, our goal will be to compute rg satisfying the error tolerance 

||r(l)-rQ||^ <6 (7.1) 
while minimizing the computational eff'ort 

Q 

subject to the constraints on {Pq} and {sq} given below, where W > is the work per iterative 
step which we scale to W = 1. We note that the preceding assumes that applying the iterative 
solver is the most computationally expensive operation and all iterations are equally expensive. 
We first formulate the optimization problem with constraints, and we then simplify the problem 
by observing that some of the inequality constraints can be replaced by equality constraints. In 
this section, we now consider a continuous, decreasing contraction radius 6{s) and a continuous, 
positive bound k(s). The load steps taken to achieve the error goal will not be uniformly spaced 
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which will take advantage of the large initial contraction region and the fact that low error is only 
desired at the endpoint, s = 1. 

We define the error at Sg by Cg = \ \r{sq) — rq\\^ for all q = 0, . . . ,Q. Then a bound on the error 
in the initial guess for g = 1, . . . , Q can be given by 

< K{Sg) - K(Sg_l) + eg_l, 

where k{s) = fc(r)dr. If ^{sg) — + e^-i < S{sg), the mapping Tg^ is a contraction and 

the error satisfies the bound 

r{sg)-4'' <a^'||r(.,)-rO||^ 

since r{sq) is a fixed point of Tg^. We assume that cq < 5(0), and we let {7q}^^Q be the supersolution 
for {e^l^^Q defined by the recurrence 

7g = a^i {K{sg) - K{sq-i) + 7g-i) , for g = 1, . . . , Q, 

70 = Co- 
in the following, we satisfy the error goal ()7.ip by making sure that the supersolution satisfies 
7Q < e- 

We now consider the question of how to achieve the error goal for the supersolution while using 
the fewest possible applications of the iterative solver. We define the set of admissible loading paths 
that satisfy the preceding theory by 

00 

= U { ({^J?=i' {^J?=o) C ' X Z>o X [0, : 
Q=l 

= So < si < ■ ■ ■ < SQ = 1, (7.2) 

K{Sg) - K(Sg_l) + 7g_l < 6 (s g) 

for all q = 1, . . . ,Q, and 7q < e|. 
Figure [6] shows the error ||r(s) — f(s)||Q^ for hypothetical loading paths 

f (S) = rg_l for S G Sg), 

using the worst-case error bound 

r{sg) -rq = a^''{r{sg) - Vg^i). 
We will next consider minimizing the work with respect to all admissible paths. 

Problem 7.1. Given e > 0, k{s) > continuous and increasing, 5{s) > continuous and decreas- 
ing, < 70 < (5(0), and < a < 1, find 

Q 

argmin y^{{Pq},{sg}) = argmin /.Pg- 

{{P<i},{s,})cA {{P,},{s,})cA 

We can see that A is non-empty by considering paths with sufficiently many small steps so that 
the error stays within the contraction domain of the iteration. Thus, the problem has a minimizer 
since the work for each path is a positive integer. We denote the minimum possible work by Wmm- 
The above problem can be analytically solved by the using following two lemmas which characterize 
paths that are optimal in the sense of this problem. We first observe that it is optimal to only 
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Figure 6. Error e{s) = | |r(s) — r(s)| |^ for the deformation problem given in 
Sectional where f(s) = r^-i, s G [sq-i,Sq). Two example admissible loading paths 
(I7.2P are displayed in bold. The contraction radius S{s) corresponding to a = | 
bounds the estimated error curve for any admissible loading path, (a) Generic path 
with multiple iterations per step, (b) Path with a single iteration per step and error 
estimate just less than d{s). Path is optimal solution for Problem 17.11 



do enough work to stay within the contraction bounds, that is, the minimum work lies on the 
boundary K{sq) - + 7g„i = 6{sq). 

Lemma 7.1. Let C = {{Pq}'^^i, {sq}^^^) € A denote an admissible loading path. Then there is 
£ = {{Pq}'^^-^, {sq}q=o) & A and J, 1 < J < Q , such that Pq = Pq for all q = I, . . . ,Q] 

K(Sq) - K{Sq-l) + %-l = 5(Sq) (7.3) 

for every q = 1, . . . , J — 1; and Sg = 1 for every q > J. Furthermore, 7q < 7q with equality if and 
only if C = C. 

The full proof is given in the Appendix. The idea is that since our goal is to only get accuracy at 
sq (|7.ip . reducing error early results in extra total work. If (j7.3p is not satisfied at some Sq, then 
we can take a larger step between Sq-i and Sq and smaller steps later (for g = 1, . . . , Q — 1) which 
reduces the supersolution for the error for all subsequent steps. 

Next, we denote the set of all admissible loading paths with work X^^i -fg ^ ^ by 

Q 

Am = { {{P,}%i, {s,}%o) ^A:Y,Pq<m] for m > 0. 

9=1 

If the minimum total work is denoted by Wmmj then Ayy^i^ is non-empty. We now show that it is 
optimal to only do one iteration per step, increasing the number of steps as necessary. 

Lemma 7.2. There is a path ({Pg}^]^, {sg}^Q) € -^Wm,i„ such that Pq = 1 and (j7.3p hold for 
all q = 1, . . . ,Q — 1. This uniquely determines Q. Furthermore, this path has the lowest estimated 
error, jq, of all paths in Aw^^^. 
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The full proof is given in the Appendix where it is shown that any step other than the last with 
Pq > 1 can be split to create a new admissible path that has the same total work and a lower error. 

Combining these two lemmas, we can characterize the optimal path for solving Problem 17.11 by 
the following algorithm, where some equations are given in implicit form: 

let 70 = eo, so = 0, q = 
while Sq < 1 

q = q+1 

Sq = SOlve(K;(Sg) - K(Sg_l) + 75_l = 6{Sq)) 



1q 



end 

Po 



mm{sq, 1) 
: 1 
aS{sq) 

loge-log(K(l)-K(sQ,i)+7Q_i) 
log a 



where [x] is the least integer not less than x. 

Figure [6}3 depicts the loading curve and optimal loading path for our example, where we directly 
use 



k{s) :-- 



r'(r)| 



dr 



in Problem 17.11 We note that the solution depicted in Figure Eb uses information about the exact 
solution, both in the growth estimate k and in the computation of contraction regions ()4.3p . In 
practice, the results will be applied using estimates to determine 6 and k. The lemmas provide the 
general intuition that it is efficient to do many steps with a single iteration per step, rather than 
fewer steps with many iterations per step. 



Appendix A. Proofs of Lemma 17.11 and Lemma 17.21 

We present detailed proofs of Lemma 17.11 and Lemma 17.2^ which use very similar estimates to 
show that a given new path is computationally favorable. 

Proof of Lemma \7.1\ Let J be the smallest integer such that sj = 1. If J = 1, then we are done; 
otherwise we show that if (|7.3p does not hold for some g = 1, . . . , J — 1, then we can adjust {sq} 
to satisfy (j7.3p with a strict decrease in the total error. 

Let j be the smallest integer such that K{sj) — K{sj-i) + 7j_i < 6{sj). If j < J, then our step 
was too conservative so we define a new stepping path {sg}^=o- > be chosen so that 

n{sj + As) — K{sj-i) + 7j_i = 6{sj). We let As = min(As, 1 — Sj) and define Sq by 

Sq, q < j, 

Sj + As, q = j, 

max(sg, Sj + As), q > j. 

This gives a new loading path = sq < si < • • • < sg = 1. By construction, steps q = 1, . . . , j 
satisfy ()7.3p . but we must show that all subsequent steps are inside the contraction region, that is 
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K{sq) — K{sg-i) + 7g_i < 5{sg) fov q > j. If > Sj + As, we have 
7i+i = a^'+'Ksj+i) - K{sj)+jj] 

= a^^+i [k(sj+i) - K{sj + As) + a^' (k(sj + As) - k(sj_i) + ij-i)] 



a 



Hk(s 



j+iJ 



K{sj + As) + a^{K{sj + As) - K{sj)) 



+ Q^^' (K(Sj) - k(Sj_i) + 7j-l)] 

= 7i+i- 

The above shows that the error supersolution 7j+i is reduced and, by consideration of the terms in 
brackets, that the approximation is inside the contraction region at Sj+i. A similar argument holds 
for the first non-degenerate step, sj* > Sj, in the case Sj + As > Sj+i. Since the remaining path is 
unchanged, we have 7^ < 7^ for all g = j, . . . , Q. We continue this process until the hypotheses are 
satisfied. □ 

Proof of Lemma \ 7.2[ We choose a path in Ayy^^^ of the form given by Lemma I7.1i Now, if 
sj = 1 for J < Q, we can combine step J and J + 1 by letting Pj = Pj+i + Pj, thus we can 
consider paths where (17. 3p holds for all q = 1, . . . ,Q — 1. 

Now, we show that if Pj > 1 for some j < Q — 1, then the error can be reduced by adding a new 
load step between Sj and s^+i. Suppose the path satisfies (j7.3p for all (7 = 1, . . . , Q — 1 and Pj > 1 
for some j < Q — 1. We will consider the new path {{Pq}'^^^, {sq}'^^^) E Awmm given by 







1 


Q = j 


Pj-^ 


q = j + l 




q> j + 1 


Sq 


q<j + l 


Sj + As 


q = j + l 


Sq-l 


q>j + l 



where As is chosen such that 

K{Sj + l) - K{Sj) +7j 

= K{sj + As) — K{sj) + a6{sj) 
= 6{sj + As) 
= ^(sj+i)- 

The above has a solution, with < As < Sj^i — Sj, by the Intermediate Value Theorem. The 
above path clearly has the same total work as the original, and we now show that it satisfies the 
contraction region constraints in Problem 17.11 We find that 

K{Sj+2) - k(Sj + i) +7j + l 



M(s 



= K{sj+2) - ^(sj+i) + a' 
= K(sj+i) — K{sj + As) 

+ Q^J"^ {K{sj + As) - K{sj) + a6{sj)) 
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Thus, we have lowered the supersolution for the error. We can apply Lemma 17.11 and the above 
argument until the hypotheses are satisfied, and at each step the error supersolution is reduced. □ 
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